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ABSTRACT 


A  graphical  test  bed  in  which  the  results  of  a  simula- 
tion experiment  can  be  reported  and  analyzed  is  described. 
The  test  bed  is  based  on  the  regression  adjusted  graphics 
and  estimation  methodology  developed  by  Heidelberger  and 
Lewis  [Ref.  1J  for  regenerative  simulation.   From  the 
graphics  and  associated  numerics,  the  experimenter  can  sum- 
marize and  see  simultaneously  relative  properties,  such  as 
bias,  normality  and  standard  deviation,  of  several  estimators 
of  a  characteristic  of  a  population  for  up  to  8  sample  sizes. 
The  evolution  of  these  properties  with  sample  size  is  also 
displayed.   The  graphics  is  supported  on  a  line  printer  to 
make  it  and  the  program  portable.   The  technique  is  illus- 
trated by  examples  concerning  the  effects  of  changes  in 
data  distribution  on  the  behavior  of  the  lag  one  serial 
correlation  coefficient,  the  estimation  of  the  shape 
parameter  of  Gamma  random  variables  and  a  comparison  of 
different  methods  (jackknife,  bootstrap)  for  estimating  the 
standard  error  of  an  estimator. 
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I.   SYNOPSIS 


SIMTBED 


THE  PROGRAM: 

Portable  FORTRAN  program  using  printer  plot  graphics 

(3  different  program  versions) 

Program  will  run  on: 

IBM 

VAX 

IBM  PC 
etc . 

ca.  9  00  lines  of  FORTRAN  Code 

Memory  requirements : 

SIMTB1      1  M  Bytes 

SIMTB2      1  M  Bytes 

SIMTB3    0 .5  M  Bytes 

(may  slightly  differ  with  different  type  of  estimator 

functions  and  subsample  sizes) 

PURPOSES: 

To  explore  the  distribution  of  a  statistical  estimator 
To  see  how  that  distribution  changes  with  sample  size 
To  compare  that  distribution  with  the  distribution  of 
competing  estimators 
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THE  USER  SUPPLIES: 

A.  Optional  Parameters: 

NE(1) ,NE(2) , . . . ,NE ( 8)   =  Subsample  Sizes  (maximum  is  8) 

The  estimator  is  computed  based  on  NE(i)  data  points 

N     =  Total  Number  of  simulated  data  points  per 
replication 
At  Subsample  Size  NE(i),  there  are  (N/NE(i)J 
independent  values  of  the  estimator 

M     =  Number  of  Replications 

When  all  replications  have  been  run,  there  are 
M*JN/NE(i)J   independent  values  of  the  estimator 
at  each  NE(i) 

D     =  Degree  of  Regression  (maximum  D  =  3) 

L     =  Number  of  Subsample  Sizes  (maximum  L  =  8) 

GRAPHICS  and  SCALING  options 

B.  Data: 

A  total  of  M*N  simulated  data  values  are  needed.   In 
SIMTB1  and  SIMTB2  the  same  data  is  used  at  each 
subsample  size  (NE(i)) .   In  SIMTB3  new  data  is  always 
generated. 

C.  Estimators: 

Up  to  3  FORTRAN  functions  (i.e.  Estimators)  are  needed. 
They  must  accept  as  inputs  a  data  subsample  and  the  size 
of  that  subsample.   They  must  return  one  value  of  the 
estimator  for  that  subsample. 
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SIMTBED  PRODUCES: 

A  one  page  graphical  output  (box  plots)  at  each 
subsample  size 

Numerical  Summaries  at  each  subsample  size  (mean, 
Std.  Dev.  ,  Std.  Dev.  of  the  mean,  skewness,  kurtosis) 
Regressions  to  quantify  changes  in  the  mean  and 
variance  as  subsample  size  changes 
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II.   INTRODUCTION 

SIMTBED,  with  the  different  versions  (SIMTB1,  SIMTB2 
and  SIMTB3)  is  a  graphical  display  program.   The  program 
is  based  on  the  program  RAGE  [Ref .  2] .   It  is  used,  with 
simulated  data,  on  a  digital  computer.   The  program  can  be 
used  to  examine  statistical  estimators  of  different  type, 
or  properties  of  a  single  estimator  under  different  distri- 
butional assumptions.   The  distribution  of  the  estimator  can 
be  explored  for  given  sample  sizes  and  the  properties  can  be 
compared  for  different  sample  sizes.   The  estimation  condi- 
tions are  controlled  by  the  experimenter.   It  is  also  possible 
to  examine  the  effects  of  changes  in  the  underlying  distri- 
bution of  the  data. 

When  the  program  SIMTBi  or  SIMTB2  is  used  with  simulated 
data,  the  data  is  assumed  to  be  independent  and  identically 
distributed  (iid).   This  iid  data  can  be  sectioned  into 
M  independent  blocks  of  specified  sample  size  N.   The  sample 
of  size  N,  will  be  sectioned  into  smaller  subsample  of  size 
NE(k).   The  estimates  are  then  calculated  from  this  subsample 
of  size  NE (k) . 

One  salient  feature  of  the  program  versions  SIMTBI  and 
SIMTB2  is  that  they  use  the  same  batch  of  simulated  random 
variables  to  explore  the  properties  of  all  the  estimators 
at  the  various  subsample  sizes.   This  is  done  for  economy 
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and  could  be  important  on  slow  computers;  the  price  paid  is 
that  the  regression  analysis  provided  by  SIMTB1  and  SIMTB2 
of  its  graphical  output  is  performed  on  correlated  samples. 

The  version  SIMTB3  uses  one  dimensional  data  and  does 
not  have  this  repetition  feature.   New  data  is  used  for  each 
calculation  of  each  estimator  at  each  subsample  size.   The 
data  is  generated  when  the  estimator  function  is  called  and 
only  the  needed  batch  of  the  exact  subsample  size  is  gener- 
ated.  This  technique  reduces  the  memory  requirements. 

Moreover  all  data  sets  are  uncorrelated  and  a  much  more 
precise  correlation  can  be  performed  if  required.  But  this 
technique  increases  the  computer  time. 

To  use  the  program  it  is  necessary  only  to  define  the 
optional  parameters  (see  Section  IV),  supply  the  simulated 
random  variables  or  a  batch  of  data  points,  and  provide  the 
FORTRAN  functions  for  the  calculation  of  the  estimators 
which  are  of  interest.   The  program  itself  will  subdivide 
the  input  data  and  feed  the  data  properly  into  the  functions, 
scale  the  graphic  display,  produce  boxplots  and  summary 
statistics.   A  regression  will  be  computed  for  the  mean  and 
variance  of  each  estimator  based  on  inverse  subsample  size. 
The  result  of  the  regression  is  displayed  graphically  and 
numerically . 

The  program  is  written  in  ANSI  Standard  FORTRAN  (X3.9- 
1966)  and  extensively  tested  on  an  IBM  3033  computer  using 
FORTRAN  IV  (H  Extended)  or  FORTRAN  IV  (Gl)  compilers.   The 
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program  SIMTBED  (all  versions)  provides  all  subroutines 
used  inside  the  program.   Besides  the  estimator  functions, 
the  user  has  to  provide  NO  additional  subroutines  via 
packages  like  IMSL.   The  program  should  be  totally  portable; 
it  has  been  tested  under  FORTRAN  77  on  a  VAX  11/780.   The 
only  limitation  is  the  available  memory. 
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III.   GENERAL  IDEA  AND  DATA  STRUCTURE  OF  THE  PROGRAM 

The  main  purpose  of  SIMTBED,  with  the  different  versions, 
is  to  explore  the  distributional  behavior  of  estimators 
and  show  their  properties  in  a  graphical  and  numerical 
display.   All  versions  use  the  same  ideas,  they  only  differ 
in  the  type  of  data  and  the  way  the  data  is  used  inside  the 
program. 
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Figure  1.   Sectioning  of  the  Data  into  M*N  Sections 
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To  study  the  behavior  of  an  estimator  the  experimenter 
usually  uses  well  known  simulated  data.   Thus  if  one  is 
interested  in  exploring  the  behavior  of  estimates  of  the 
shape  parameter  in  a  Gamma  population,  one  generates  Gamma 
variates  from  a  random  number  generator  package  (e.g.,  IMSL 
subroutine  Chap.  G) .   This  batch  of  simulated  data  is 
processed  by  SIMTBED  in  the  following  way. 

The  data  batch  is  first  divided  into  M  independent 
blocks.   Each  block  contains  N  data  points.   So  the  starting 
data  batch  has  to  consist  of  M*N  data  points  (see  Figure  1) . 

All  blocks  are  divided  into  subsamples  of  size  n..   The 

r  1 

actual  subsample  size  n.  is  an  element  of  the  subsample 
size  array  NE .   This  array  can  store  up  to  8  different 
values .   Then  the  estimator  is  calculated  for  each  subsample 

of  size  n. .   The  estimator  function  will  be  calculated 

i 

( [N/n . J ) *M  times.   This  total  population  of  estimates  is  used 
to  evaluate  the  summary  statistics  for  the  estimate  and 
construct  the  corresponding  box  plot.   If  NE  contains 
another  element,  the  blocks  are  divided  into  the  new  sub- 
sample  n.  ,  size  and  all  calculations  are  done  again. 
In  addition  to  the  summary  statistics  and  the  box 
plots  (see  e.g.,  Figure  3a),  a  regression  on  the  averages 
and  on  the  variance  is  computed,  following  the  methodology 
of  regression  adjusted  estimate  (RARE)  developed  by 
Heidelberger  and  Lewis  [Ref .  1] . 
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The  RARE  estimate  is  the  regression  coefficient  an . 
It  is  the  asymptotic  estimate  of  the  expected  value  of  the 
parameter.   The  unbiased  RAGE  estimate  of  the  average  of  the 
parameter  is  determined  by  the  regression  formula: 

E<e(n.))   =   aQ  +  a±   i  +  «2  \   +  ...  +  «D  \ 

1       n .  n . 

1  1 

where  D  is  an  input  parameter. 

The  asymptotic  RARE  estimate  is  printed  as  a  dashed  line 
in  the  graph. 

The  unbiased  RAGE  estimate  of  the  variance  of  the  esti- 
mator is  given  by  the  formula: 


S2(n  )   =   a   -L  +  0   L_  +  6_  — i-y  + 

l        In.     2/\-L.j     3/\^ 
i       (n. )        J  (n. ) 


The  regression  calculations  of  the  average  are  done  for 
each  block  separately.   The  result  are  M  sets  of  regression 
coefficients.   The  finally  printed  regression  coefficients 
are  the  averages  of  the  M  replications  (see  Figure  2) .   From 
the  set  of  coefficients  the  variance  and  the  standard 
deviation  of  the  regression  coefficients  is  calculated. 

The  regression  on  the  variance  is  done  once,  using  all 
(M*N)  data  points. 
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Figure  2.   Structure  of  the  Regression  Coefficients 
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IV.   ARGUMENTS  OF  THE  PROGRAMS  SIMTB1  AND  SIMTB2 

All  versions  of  SIMTBED  have  in  general  the  same  argu- 
ment list.   SIMTB2  (a  version  for  multivariate  random  varia- 
bles) has  two  additional  arguments.   All  arguments  will  be 
described  in  detail  and  all  restrictions  or  limitations  will 
be  discussed.   In  Section  V,  detailed  examples  of  setups 
for  SIMTB1  and  SIMTB2  are  given. 

A.  X — DATA  ARRAY 

The  array  X  is  the  input  data  array.   X  is  single  pre- 
cision REAL*4,  and  is  generally  simulated  data  e.g.,  in 
Section  VII,  independent  Gamma  variates. 

In  SIMTB1  the  dimension  of  the  array  X  must  not  exceed 
50,000. 

In  SIMTB2  the  input  data  is  multivariate  and  conse- 
quently the  array  X  has  two  dimensions.   The  size  of  the 
dimensions  is  not  directly  limited  by  the  program.   The 
space  must  be  provided  by  the  calling  program  and  is  passed 
as  an  argument  (see  IR  and  IRK) .   The  memory  requirements 
increase  rapidly  as  the  dimensions  of  X  increase. 

B.  N — SAMPLE  SIZE 

The  sample  size  N  is  the  number  of  data  points  per  section 
of  input  data  X.   N  is  an  INTEGER.   Depending  on  the  pre- 
cision of  the  simulation  that  is  required,  the  sample  size 
N  can  vary  from  1  to  50,000. 
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In  SIMTB2  N  is  the  number  of  multivariate  data  points  per 
section.   M  (the  number  of  replications)  times  N  must  not 
exceed  IR,  the  row  dimension  of  X. 

If  M  times  N  exceeds  50,000  an  error  message  will  be 
printed  and  the  execution  will  be  terminated.   If  the  product 
M  times  N  exceeds  the  total  number  of  data  points  provided 
by  the  user,  NO  error  message  will  signal  the  user  error  and 
the  result  of  the  execution  is  not  predictable. 

C.  M— NUMBER  OF  REPLICATIONS 

The  number  of  replications  M  of  the  array  X  is  an 
INTEGER.   M  determines  the  number  of  sections,  into  which 
the  data  set  X  is  divided.   So  M*N  is  the  dimension  of  X  in 
SIMTB1.   The  parameter  M  also  determines  the  number  of 
regressions  on  the  mean  values  that  will  be  run  to  find  the 
regression  coefficients  in  the  regression  on  the  mean.   If 
M  is  1  only  one  regression  will  be  done  and  no  variance  and 
standard  deviation  of  the  regression  coefficients  can  be 
calculated. 

D.  NE — SUBSAMPLE  SIZE  ARRAY 

The  argument  NE  is  an  INTEGER  array,  containing  up  to 
8  subsample  sizes.   These  are  the  subsample  sizes  at  which 
the  properties  of  the  estimator  are  to  be  investigated.   NE 
must  always  contain  8  elements.   If  less  than  8  subsample 
sizes  are  used,  the  array  must  be  filled  up  with  dummy 
values.   The  estimator  is  calculated  for  each  subsample 
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individually.   The  used  values  of  NE  (not  the  dummy  values) 
have  to  be  ordered  from  the  smallest  to  the  largest 
(NE(1)  <  NE(2)  <  NE(3)  <  . . .) . 

The  program  can  only  handle  up  to  a  total  number  of 
12,500  estimates.   This  limit  is  caused  by  storage  limita- 
tions and  can  be  expanded  for  larger  computers.   The  smallest 
subsample  size  NE(1)  produces  the  most  estimates  and  if  this 
case  exceeds  12,500  the  execution  will  be  terminated  and  an 
error  message  will  be  printed. 

In  SIMTB2  all  values  of  NE  have  to  be  less  than  5,000. 

E.  L— NUMBER  OF  SUBSAMPLE  SIZES  (BOX  PLOTS) 

The  number  of  box  plots  L  in  the  graphical  output  (the 
number  of  subsamples  at  which  the  estimator  is  examined) , 
is  an  INTEGER  with  values  from  1  to  8 .   The  value  of  L 
determines  how  many  elements  of  the  array  NE  the  program 
will  use  (e.g.,  L  =  2,  the  program  uses  only  NE(1)  and  NE(2) 
for  the  calculations) .   L  determines  the  number  of  box  plots 
and  the  number  of  corresponding  summary  statistics  in  the 
output.   If  L  is  out  of  range  the  program  will  terminate 
execution  and  print  an  error  message. 

For  L  =  1  no  regression  is  possible.   No  regression 
output  will  be  printed. 

F.  D — DEGREE  OF  REGRESSION 

The  degree  of  regression  on  the  mean  D,  is  an  INTEGER 
with  values  from  1  to  3 .   The  chosen  degree  refers  to  the 
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number  of  coefficients  calculated  and  printed  for  the 
regression  equations.   Experience  has  shown  that  D  =  3  is 
preferable  (higher  values  cause  severe  numerical  problems  in 
the  regression  computations). 

If  D  exceeds  the  value  of  L-l,  the  program  reduces  D 
to  this  number,  regardless  of  the  value  chosen  by  the  user. 

G.   RG — REDUCED  GRAPHICS 

The  argument  RG  is  an  INTEGER  with  value: 

0  =  off  ==>  full  graphics 

1  =  set  ==>  reduced  graphics 

In  the  reduced  graphics  case  the  vertical  scale  of  the 
graph  is  reduced  by  eliminating  the  outliers.   Only  the 
number  of  outliers  is  counted  and  the  number  is  printed. 
This  enables  the  user  to  see  the  body  of  the  boxplot  in  more 
detail . 

H.   SEI — SCALING  ESTIMATORS  INDIVIDUALLY 

The  argument  SEI  is  an  INTEGER  with  value: 

0  =  off  ==>  common  scaling 

1  =  set  ==>  individual  scaling. 

For  SEI  =  0  all  printed  graphs  (max.  3  per  program 
run)  are  scaled  to  the  same  range.   This  makes  the  compari- 
son of  the  different  graphs  easier.   The  value  SEI  =  1  causes 
the  program  to  scale  each  graph  individually. 

The  argument  SEI  works  in  combination  with  the  arguments 
RG  and  SVS .   The  combination  of  these  three  arguments  make  it 
possible  to  fit  the  printed  graphs  to  the  needs  of  the  user. 
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I.   SVS--SETTING  THE  VERTICAL  SCALE 

The  argument  SVS  is  an  INTEGER  with  value: 

0  =  off  ==>  automatic  scaling 

1  =  set  ==>  scaling  by  the  user. 

SVS  =  0  causes  the  program  to  calculate  the  scale  of 
the  graphic  printout.   The  final  graphic  display  is  influ- 
enced by  the  chosen  values  of  RG  and  SEI .   The  values  of  YMIN 
and  YMAX  are  ignored  by  the  program. 

For  SVS  =  1  the  user  must  provide  the  scaling.   The 
graphics  are  vertically  scaled  to  a  given  minimum  value 
(YMIN)  and  maximum  value  (YMIN  and  YMAX) .   The  arguments  of 
RG  and  SEI  are  ignored.   Only  the  numbers  of  outliers  outside 
the  display  are  printed  for  each  box  plot. 

J.   YMIN — MINIMUM  VALUE  OF  THE  VERTICAL  SCALE 

The  scaling  parameter  YMIN  is  data  type  REAL*4.   It  is 
the  lower  limit  of  the  vertical  scale.   It  affects  the 
scaling  only  in  combination  with  SVS  =  1.   If  the  chosen 
value  is  so  large  that  the  value  YMIN  lies  inside  the  body 
of  the  box  plot,  an  error  occurs  and  the  program  ends  with 
an  abnormal  ending. 

K.   YMAX — MAXIMUM  VALUE  OF  THE  VERTICAL  SCALE 

The  scaling  parameter  YMAX  is  data  type  REAL*4.   It  is 
the  upper  limit  of  the  vertical  scale.   The  scaling  is  only 
effected  if  SVS  has  value  1.   If  the  value  of  YMAX  is  too 
small  and  it  lies  inside  the  body  of  the  box  plot,  the 
program  comes  to  an  abnormal  ending. 
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YMIN  and  YMAX  in  combination  with  SVS  =  1  should  only 
be  used  for  well  known  graphic  ouptut.   With  this  option, 
it  is  possible  to  scale  the  output  of  different  program 
runs  to  a  common  scale.   In  particular  if  more  than  3  esti- 
mators have  to  be  estimated  and  compared,  so  that  a  common 
scale  is  needed,  this  option  may  be  used. 

L.   NEST — NUMBER  OF  ESTIMATORS 

The  parameter  NEST  is  an  INTEGER  with  the  value  1,  2 
or  3.   The  value  of  NEST  determines  the  number  of  different 
one-page  graphic  displays  the  program  produces  with  one  call 
from  the  calling  program.   Usually  the  value  of  NEST  is 
equal  to  the  number  of  different  estimators  used  in  the 
program. 

In  SIMTB2  the  same  estimator  may  be  used  with  different 
(e.g.,  normal,  exponential  etc.)  distributed  data  sets. 

M.   EST1,  EST2,  EST3 — ESTIMATOR  FUNCTIONS 

These  3  parameters  are  data  type  REAL*4  and  are  used 
to  pass  the  EXTERNAL  function  names  to  the  program.   An 
external  declared  function  is  a  function  which  computes  the 
value  of  an  estimator.   This  function  may  call  other  routines, 
but  the  final  value  of  the  estimator  must  be  passed  over 
this  function  name. 

For  SIMTBl  each  function  must  have  the  two  arguments  X 
and  NEK  (e.g.,  FUNCTION  VARIANCE  (X,NEK) ) .   X  is  the  data 
array  and  NEK  is  the  number  of  data  elements  in  X. 
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SIMTB2  needs  for  each  function  four  arguments  X,  NEK, 
IDR,  IRK  (e.g.,  FUNCTION  CORRELATION  (X, NEK, IDR, IRK) ) .   X 
is  a  two  dimensional  data  array  and  NEK  is  the  subsample 
size,  for  which  the  function  is  evaluated.   IDR  and  IRK  are 
the  dimensions  of  the  array  X. 

If  the  user  wants  to  use  less  than  3  estimator  functions, 
he  must  use  dummy  arguments  and  choose  the  correct  value  of 
NEST.   The  easiest  way  to  do  this  is  to  use  a  function  of 
a  previous  used  estimator  again  (for  details  of  the  pro- 
gramming see  Section  V) . 

N.   TTL1,  TTL2,  TTL3- DESCRIPTION  OF  THE  ESTIMATORS 

These  arrays  are  used  to  pass  titles  from  the  calling 
program  to  SIMTB.   Each  title  has  to  be  120  characters  long 
(15  fields,  each  8  characters  wide).   The  title  is  printed 
below  the  output  of  the  corresponding  function. 

If  the  user  uses  less  than  3  functions  and  titles,  he 
has  to  use  dummy  arguments  (for  details  see  Section  V) . 

If  the  titles  are  not  in  the  correct  format,  correspond- 
ing to  the  FORMAT  statement  the  program  will  not  execute 
properly . 

0.   IR,  IRK — DIMENSIONS  OF  THE  ARRAY  X 

These  arguments  are  used  ONLY  in  SIMTB2.   IR  Crow  dimen- 
sion) and  IRK  (column  dimension)  are  INTEGERS  and  have  to 
match  the  dimension  of  X  in  the  calling  program. 

If  the  dimensions  don't  match,  NO  error  message  will  be 
produced  and  the  output  is  unpredictable.   The  error  may 
not  be  obvious,  and  the  ouput  may  look  reasonable. 
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V.   DUMMY  EXAMPLES  OF  IMPLEMENTING  SIMTB  AND  SIMTB2 

INTO  A  DRIVER-PROGRAM 

The  following  examples  show  how  the  programs  SIMTB1  or 
SIMTB2  can  be  implemented  into  a  FORTRAN  driver-program. 
The  only  purpose  of  the  examples  is  to  clarify  the  FORTRAN 
implementation,  to  avoid  programming  errors  by  the  user. 
Additional  comment  lines  are  added  to  the  program  examples 

In  the  first  example  SIMTB1  is  used  to  compare  two 
different  estimators  (VAR  and  UNBVAR)  of  the  variance  of  a 
normal  sample.   The  sample  may  be  generated  with  a  random 
number  generator  (e.g.,  LLRANDOMII) . 

In  the  second  example  SIMTB2  is  used  to  compare  the 
distribution  of  two  estimators.   The  estimators  are  the 
covariance  (COV)  and  the  correlation  coefficient  (CORR)  of 
a  bivariate  standard  normal  sample. 

A.   EXAMPLE  1  USING  SIMTBl 

MAIN 
C   EXAMPLE  of  SIMTBl    Calling  program,  it  has  not  to  be 
the  MAIN  program 

REAL*4  X(50000),  YMIN,  YMAX,  VAR,  UNBVAR 
REAL*8  Tl(15),  T2(15),  T3(15) 

INTEGER  N,  M,  NE(8),  L,  D,  RG,  SEI,  SVS ,  NEST 
C 

DATA  N    /2  50  0/ 
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DATA  M    /   20/ 

DATA  NE   710,20,25,50,100,0,0,0/ 
C  Array  NE  must  have  8  elements,  if  only  5  are  used,  add 
C   dummy  variables  and  set  L  =  5 

DATA  L    /    5/ 

DATA  RG   /    0/ 

DATA  SEI  /    0/ 

DATA  SVS  /    0/ 

DATA  NEST/    2/ 

DATA  Tl   /'ESTIMATE1 , 'OF  THE  V ', 'ARIANCE ',' USING   ', 
+'VAR=(1/N) ' , '*SUM(X(I' , * )-XBAR) *' , ' *2       ',7*'  '/ 

DATA  T2  /'ESTIMATE' , 'OF  THE  V',' ARIANCE  ',' USING   ', 
+'VAR=(1/(1' , '-N) ) *SUM' , ' (X(I)-XB' , 'AR) **2   ',7*'  '/ 
C  All  15  fields  (each  8  characters)  have  to  be 
C   initialized 
C 

EXTERNAL  VAR,  UNBVAR 
C 

C   Generate  M*N  independent  Normal  (0,1)  Random  numbers 
C  and  store  into  X 
C 

CALL  S IMTB1 ( X , N ,M , NE , L , D , RG , SEI , SVS , YMIN , YMAX , NEST 
,VAR,Tl,UNBVAR,T2,VAR,Tl) 
C   EST3=VAR  and  TTL3=T1  used  as  dummy  variables 
C 

STOP 

END 
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c 

FUNCTION  VAR  (X,NEK) 
C   Function  to  calculate  the  Variance. 

C  All  calculations  inside  the  function  should  be  done  in 
C   DOUBLE  PRECISION 
C 

REAL*4  X(N) ,  VAR 
REAL* 8  SUM,  XBAR ,  DVAR 
C 

SUM=0.0D0 
DO  10  1=1, N 

SUM=SUM+DBLE(X(I) ) 
10  CONTINUE 

XBAR=SUM/FLOAT (N) 
SUM=0.0D0 
DO  20  1=1, N 

SUM=SUM+( (DBLE(X(I) ) ) -XBAR) **2 
20  CONTINUE 

DVAR= SUM/ FLOAT (N) 
VAR=SNGL(DVAR) 
C 

RETURN 
END 

FUNCTION  UNBVAR  (X,NEK) 
C   Function  to  calculate  the  Variance. 
C  All  calculations  inside  the  function  should  be  done  in 
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C   DOUBLE  PRECISION 
C 

REAL*4  X(N),  UNBVAR 
REAL* 8  SUM,  XBAR,  DUNVAR 
C 

SUM=0.0D0 
DO  10  1=1, N 

SUM=SUM+DBLE(X(I) ) 
10   CONTINUE 

XBAR= SUM/FLOAT (N) 

SUM=0.0D0 

DO  20  1=1, N 

SUM=SUM+( (DBLE(X(I) ) ) -XBAR) **2 
20   CONTINUE 

DUNVAR= S UM/FLOAT ( N ) 
UNBVAR=SNGL ( DUNVAR) 
C 

RETURN 
END 

B.   EXAMPLE  2  USING  SIMTB2 

MAIN 
C  EXAMPLE  of  SIMTB2    Calling  program,  it  has  not  to  be 
the  MAIN  program 
REAL*4  X(25000,2),  YMIN,  YMAX,  COV,  CORR 
REAL* 8  Tl(15),  T2(15),  T3(15) 
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INTEGER  N,  M,  NE(8),  L,  D,  RG,  SEI,  SVS,  NEST,  IR, 
+  IRK 

C 

DATA  N    /250  0/ 

DATA  M    /   10/ 

DATA  IR   /25000/ 

DATA  IRK  /    2/ 
C   IR  and  IRK  must  be  equal  to  the  dimensions  of  X 

DATA  NE   710,20,25,50,83,100,125,250/ 

DATA  L    /    8/ 

DATA  RG   /    0/ 

DATA  SEI  /    0/ 

DATA  SVS  /    1/ 

DATA  YMIN/  0.0/ 

DATA  YMAX/  1.0/ 

DATA  NEST/    2/ 

DATA  Tl   /'ESTIMATE' , 'OF  THE  C ' , ' OVARIANC * , ' E 
+11*'   '/ 

DATA  T2  /'ESTIMATE' , 'OF  THE  C ' , ' ORRELATl ' , 'ON  COEFF ' , 
+'ICIENT    ' , '10*'  '/ 
C  All  15  fields  (each  8  characters)  have  to  be 
C   initialized 
C 

EXTERNAL  COV,  CORR 
C 

C   Generate  M*N  pairs  of  independent  random  bivariate 
C  numbers,  each  pair  being  independent  N(0,1)  random 
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C   variables  and  store  into  X 
C 

CALL  SIMTB2  (X ,N ,M,NE, L, D, RG , SEI , SVS , YMIN , YMAX ,NEXT 
+  ,C0V,Tl,C0RR,T2,C0V,Tl,IR,IRK) 

C   EST3=C0V  and  TTL3=T1  used  as  dummy  variables 
C 

STOP 

END 
C 

c 

FUNCTION  COV  (X, NEK, IDR, IRK) 
C    Function  to  calculate  the  Covariance. 

C   All  calculations  inside  the  function  should  be  done  in 
C    DOUBLE  PRECISION 
C 

REAL*4  X( IDR, IRK),  COV 

REAL* 8  SUM1,SUM2,SUM3,XBAR1,XBAR2,EX1X2,DC0V 
C 

SUM1=0.0D0 
SUM2=0.0D0 
SUM3=0.0D0 
DO  10  1=1, N 

SUM1=SUM1+DBLE ( X ( 1 , 1 )  ) 
SUM2=SUM2+DBLE (X ( I , 2 ) ) 
SUM3=SUM3+DBLE(X(I,1) *X(I,2) ) 
10  CONTINUE 

XBAR1=SUM1/ FLOAT (N) 
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XBAR2=SUM2/FLOAT (N) 
EX1X2=SUM3/FL0AT (N) 
DC0V=EX1X2- (XBAR1*XBAR2) 
COV=SNGL(DCOV) 
C 

RETURN 
END 

C 

FUNCTION  CORR  (X,NEK, IDR, IRK) 
C  Function  to  calculate  the  Correlation  coefficient 
C  All  calculations  inside  the  function  should  be  done  in 
C   DOUBLE  PRECISION 
C 

REAL*4  X(IDR,IRK),  CORR 

REAL*8  SUM1 , SUM2 , SUM3 , XBARl , XBAR2 , EX1X2 ,  VARl , VAR2 , 
+        COV,DCORR 
C 

SUM1=0.0D0 
SUM2=0.0D0 
SUM3=0.0D0 
DO  10  I=1,N 

SUM1=SUM1+DBLE(X(I/1)  ) 
SUM2=SUM2+DBLE(X(I,2)  ) 
SUM3=SUM3+DBLE(X(I,1) *X(I,2) ) 
10   CONTINUE 

XBARl=SUMl/FLOAT (N) 
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XBAR2-SUM2/FLOAT (N) 

EX1X2=SUM3/FL0AT (N) 

SUM1=0.0D0 

SUM2=0.0D0 

DO  20  1=1, N 

SUM1=SUM1+DBLE(X(I,1) **2) 
SUM2=SUM2+DBLE(X(I,2) **2) 
20   CONTINUE 

VARl=(SUMl/FLOAT(N) )- (XBARl**2) 

VAR2=(SUM2 /FLOAT (N) )-(XBAR2**2) 

C0V=EX1X2- (XBAR1*XBAR2) 

DCORR=COV/( (VAR1*VAR2) **0.5) 

CORR=SNGL(DCORR) 

RETURN 
END 
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VI.   STUDY  OF  THE  BEHAVIOR  OF  SERIAL  CORRELATION 
ESTIMATES  FOR  DIFFERENT  DISTRIBUTIONS 

A.  CALCULATION  OF  THE  FIRST  SERIAL  CORRELATION  COEFFICIENT 
It  is  known  that  for  an  independent  sample  from  a 

population  with  finite  variance,  the  distribution  of  the 
serial  correlation  coefficient  (Anderson  and  Walker,  1964) 
[Ref.  3]  is  asymptotically  Normal  with  mean  zero  and  vari- 
ances 1/n,  where  n  is  the  sample  size.   If  the  population  is 
i.i.d  Normal  then  the  bias  is  exactly  -1/n.   Since  those 
asymptotic  properties  are  frequently  used  as  approximations 
in  tests  of  significance,  it  is  important  to  know  how  valid 
the  approximation  would  be  in  small  samples  from  a  variety 
of  distributions.   We  will  look  at  that  question  in  the  next 
two  sections  and  then  go  on  to  consider  two  alternative 
measures  of  correlation,  Fisher's  z-transform  and  the  2-fold 
jackknifed  estimate  of  the  correlation.   Their  ability  to 
reduce  bias  and/or  induce  Normality  will  be  examined  against 
other  changes  in  the  distribution  of  the  estimators,  particu- 
larly variance  inflation.   A  simulation  study,  without 
graphics,  of  some  of  these  problems  was  conducted  by  Cox 
(1966)  [Ref.  4] . 

B.  SIMTB1  OUTPUT  FOR  SERIAL  CORRELATION 

Figure  3(a)  shows  the  simulated  distribution  and  sample 
properties  of  the  serial  correlation  coefficient  estimate 


36 


r 


n  *j         (X  -Xl)(X    -Xn) 

=    i=1       J : 

n  n  - 

(n-1)  I       (X   -  X  )Z 
j-1    3 


where: 


n 


X    =    J   x./n  , 

j=l   3 
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1      j  =  l    3 


Xn   =    E   X  /(n-1) 

n      j  =  2   J 


for  various  sub-sample  sizes  n  =  n..   This  definition  matches 
that  used  by  Anderson  and  Walker  (1964).   We  consider  first 
subsamples  of  size  n,  =  10,  and  then  of  size  n^  =  20,  n^  =  30, 
n.  =  40,  n_  =  50,  nfi  =  75,  n_  =  100  and  n_  =  150,  successively 
For  each  subsample  size  the  input  sample  of  N  =  5000  simulated 
Normal  (0,1)  random  variables  is  divided  into  as  many  full 
subsamples  of  size  n.  as  possible,  and  the  serial  correlation 
is  computed  for  each  of  the  iN/n . J  subsamples  of  size  n.. 
The  entire  procedure  is  then  replicated  M  =  10  times,  each 
time  with  a  new  simulated  sample  of  N  =  5000  Normal  (0,1) 
variables . 

After  all  M  replications  have  been  run,  all  the  estimates 
of  serial  correlation  for  each  subsample  size  are  grouped 
together  and  their  simulated  distribution  is  presented  via  a 


37 


boxplot  and  summary  statistics.   The  boxplot  follows  the 
standards  discussed  in  Mosteller  and  Tukey  (1977)  [Ref.  5] 
with  the  median  denoted  by  a  +  within  the  box,  the  mean  by  a 
*  within  the  box,  the  outliers  by  O's,  and  the  far  outliers 
by  *'s  beyond  the  whiskers.   The  summary  statistics  include 
the  sample  mean,  sample  standard  deviation,  estimated  standard 
deviation  of  the  sample  mean  (i.e.,  sample  standard  deviation/ 
sqrt (M  |N/n. J  )  ,  sample  skewness  and  sample  kurtosis  of  the 
correlation  estimates. 

Looking  at  the  output,  the  first  (leftmost)  boxplot  in 
the  graph  in  Figure  3(a)  shows  the  distribution  of 


(#  Replications 


, Simulation  , 
Sample  Size 


, Subsample, 
[        Size  ' 


=   10  x 


5000 
10 


=   10  x  500  =  5000 


estimates  of  serial  correlation  from  independent  subsamples 
of  size  n,  =  10.   Summary  statistics  for  the  boxplot  can  be 
found  below  the  graph  in  the  column  labeled  "Subsample  Size 
10,"  so  that  the  average  serial  correlation  is  -.1074,  and 
the  estimated  standard  deviation  is  .299  6.   The  estimated 
standard  deviation  of  the  serial  correlation  estimate  is 


.2996//(5000)  =  .00424.   Recall  that  this  refers  to  correla- 
tion estimates  based  on  subsamples  of  size  10. 

Since  the  X-axis  of  the  graph  represents  subsample  size, 
the  last  (rightmost)  boxplot  shows  the  distribution  of 
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10  x 


5000 


150 


=   10  x  33   =   330 


estimates  of  serial  correlation  from  independent  subsamples 
of  size  n8  =  150.   Although  the  330  estimates  are  independent 
of  each  other ,  they  are  not  independent  of  the  5000  estimates 
that  comprise  the  first  boxplot  since  the  same  data  (divided 
and  processed  in  different  ways)  was  used  for  both.   Summary 
statistics  show  that  the  average  correlation  has  dropped  to 
-.007372,  indicating  the  fall  off  in  bias,  and  the  standard 
deviation  has  dropped  to  .07822,  indicating  the  greater 
accuracy  with  which  the  correlation  can  be  estimated  when  150 
points,  rather  than. 10,  are  available. 

In  order  to  quantify  the  changes  that  are  occurring  in 
the  mean  and  variance  of  the  distribution  of  the  estimator 
as  subsample  size  changes,  SIMTBl  performs  two  types  of 
regressions.   The  first  regression  is  on  the  averages  and  is 
done  after  each  replication,  using  the  average  serial  corre- 
lation for  that  replication,  r   ,  as  the  dependent  variable. 

i 
Inverse  powers  of  the  subsample  size  serve  as  the  indepen- 
dent variables.   For  Figure  3(a)  the  degree  of  the  regression 
was  chosen  to  be  D  =  3  so,  for  each  replication,  the  equations 
we  attempt  to  fit  by  least  squares  are: 


-  al    a2    a3 

rn.   =   a0  +  57  +  T  +  —    for    i  =  1,2,...,8. 
i  i    n.    n . 

l     i 

This  form  anticipates  the  general  asymptotic  expansion 


39 


a,    a«  a, 

E(0(n)  )   =  Q   +   -L+-±  +  -±+.,. 

n     2     3 

n  n 


which  holds  true  in  the  current  situation  with  0=0  and 
(in  the  Normal  case)  a  =  -1  (see  Cramer  (1948)  for  general 
results  of  this  type)  [Ref.  6]. 

Values  of  a  ,  a,,  a„ ,  and  a~  are  calculated  after  each 
replication,  averaged  across  the  M  replications  to  get  a., 
a,/  a»,  and  a  ,  and  then  the  averages  are  reported  below  the 
summary  statistics  on  the  line  "Mean  of  Regression  on  Averages- 
Coefficients."   We  find  that  a   =  -.000272  and  a   =  -1.03074, 
both  close  to  their  theoretical  counterparts. 

Because  we  have  10  replications  and  therefore  10  inde- 
pendent values  of  each  of  a~  ,  a.  ,  a„  ,  and  a..,  we  can  also 
estimate  the  variances  and  standard  deviations  of  a_,  a,,  a„, 
and  a~  across  replications.   These  values  are  presented  on 
the  two  lines  immediately  below  the  coefficients.   For 
instance,  the  estimated  s.d.  of  the  estimate  an  =  -.000272 
of  aQ  is  .003892. 

The  regression  line  for  the  mean  value  of  the  estimator 
is  presented  visually  in  the  graph  as  a  dotted  curve.   The 
estimated  asymptote  (i.e.,  a"n)  is  printed  with  a  dashed  line 
wherever  it  does  not  coincide  with  the  regression  line.   Bias, 
therefore,  can  be  viewed  as  the  difference  between  those  two 
lines. 

The  second  regression  referred  to  above  is  done  after  all 
replications  have  been  run  and  the  variances  of  the  estimators 
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at  each  subsample  size  have  been  calculated.   (Note  that  the 
standard  deviations,  not  the  variances,  are  presented  in  the 
summary  statistics.)   It  should  be  recalled  from  previous 
discussion  that  these  variances,  as  well  as  all  measures  in 
the  summary  statistics,  are  based  on  the  grouping  together 
of  the  serial  correlations  from  all  replications ,  at  each 
subsample  size.   This  is  in  contrast  to  the  procedure  for  the 
regression  on  the  means,  where  average  correlations  are  com- 
puted for  each  subsample  size  for  each  replication.   In  the 
case  of  the  variances,  we  have  8  equations: 


b0     bl  b2     b3 

Var(rn.)   =   57  +  -172  +  T  +  -572  '   i  =  l,2,...,8, 

1        1    n .  n .    n . 

1  11 


which  we  fit  by  least  squares  in  order  to  estimate  the  coeffi- 
cient 6n,  $,  ,  62/  an<^  3-}  in  tne  presumed  asymptotic  expansion 

6Q    ex    e2    63 
Var(0(n))   =  JL  +  +  -^  +  +    .  .  . 

n  '  n    n 

This  expansion  holds  for  the  variance  of  the  estimated 
serial  correlation  coefficient  for  independent  data.   Usually 
it  will  be  3Q  in  which  we  are  most  interested  since  3Q  is 
used  in  computing  asymptotic  relative  efficiencies  of  esti- 
mators.  For  independent  data  with  finite  variance,  we  know 
that  6n  =  1.   The  computed  values  of  b_,  b, ,  b2,  and  b^, 
are  presented  on  the  line  labeled  'Regression  on  Variance — 
Coefficients'.   Notice  that  bQ  =  .7438  is  close  to  the 
theoretical  value  of  1. 
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The  final  two  numbers  on  Figure  3(a),  YMIN  and  YMAX, 
simply  show  the  scale  of  the  vertical  axis.   Because  the  SIMTB1 
program  option  to  put  Figures  3(a),  3(b)  and  3(c)  on  the 
same  scale  was  in  effect,  it  may  be  that  no  boxplot  in  a 
given  Figure  (e.g.,  Figure  3(b))  requires  the  full  range  of 
Y- values . 

In  order  to  produce  Figure  3(b),  the  Normal  (0,1)  data 

that  went  into  Figure  1(a)  was  squared  to  create  longer 

2 
tailed  x  (D  random  variables.   The  output  is  entirely  analo- 
gous to  that  for  Figure  3(a) .   Similarly,  for  Figure  1(c) , 
the  Normal  (0,1)  data  was  exponentiated  in  order  to  create 
Lognormal  (0,1)  data  and  to  produce  analogous  graphical 
output.   The  indication  is  that  the  distribution  of  the  sample 
serial  correlation  is  robust  with  respect  to  the  population 
distribution. 

The  features  of  the  SIMTBl  output  will  become  clearer  when 
they  are  associated  with  the  various  properties  of  the  corre- 
lation estimator.   First,  however,  a  few  technical  comments 
concerning  the  regressions  are  necessary. 

C.   SOME  COMMENTS  ON  THE  REGRESSIONS 

Two  types  of  problems,  numerical  and  statistical,  can 
occur  when  attempting  to  fit  the  two  sets  of  regression 
equations  presented  in  Section  VI. B. 

First,  there  is  the  question  of  numerical  stability  when 

-1    -2    -  3 
the  independent  variables,  {1,  n.  ,n.  ,n.  }  or 

{n.  ,  nT  '     .    n.  ,  n~  '  }  decrease  geometrically.   If  we 
l     i    '   i  '   i 
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T 
attempt  to  form  X  X,  where  X  is  the  respective  design  matrix 

T  . 
and  X   is  the  transpose  of  X,  we  get  values  that  range  from 

8    _6 
8  (assuming  8  subsample  sizes)  to   £   n.   for  the  regression 

8    -2      8    i=1  ^ 
on  the  means,  and   J"   n.   to   7   n~5  for  the  regression  on 

1=1        1=1 

the  variances.   Experience  has  shown  that  attempts  to  solve 

T 
systems  with  such  extremes  in  the  X  X  matrix  produce  errone- 
ous results.   Instead,  SIMTB1  scales  the  design  matrices  by- 
multiplying  each  column  of  X  by  Max(n.)  raised  to  the  proper 
power  so  that  no  entry  becomes  too  small.   The  standard 
Choleski  decomposition  is  then  used  to  fit  the  equations, 
and  the  coefficients  are  properly  rescaled  before  they  are 
reported.   This  procedure  produces  numerically  reliable 
results. 

The  second  problem  concerns  the  breakdown  of  statistical 
assumptions  in  our  regression  models.   It  has  already  been 
pointed  out  in  Section  VI .B  that  the  two  sets  of  dependent 
variables : 

(1)   the  0(n.)  when  considering  the  regression  on  the  means; 


(2)   the  s  (n. )  =   L  7 

l       . L 


M  fcj/nj   (0.  (ni)  -  0(n, 


when  considering  the  regression  on  the  variances, 
are  not  independent  over  i  since  all  are  based  on  the  same 
simulated  data.   The  extent  of  the  dependence  is  demonstrated 

by  the  correlation  matrix  in  Table  1.   Entries  in  that  table 

2  2 

show  the  estimated  correlation  between  s  (n.)  and  s  (n.)  for 
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TABLE  1 

ENTRIES  IN  THE  TABLE  ARE  THE  ESTIMATED  CORRELATIONS 
BETWEEN  THE  ESTIMATED  VARIANCES  OF  THE  r    AT 
DIFFERENT  SUBSAMPLE  SIZES:  ni 

CORR(s2(rn  )/s2(rn  ))   for   i  =  1,...,8,   j  =  1,...,8 
i      nj 


1 

1.00 

.49 

.46 

-.26 

.18 

-.17 

.14 

.01 

2 

.49 

1.00 

.40 

.55 

.11 

.38 

.38 

-.03 

3 

.46 

.40 

1.00 

.23 

.23 

.44 

.21 

.29 

4 

-.26 

.55 

.23 

1.00 

.42 

.86 

.57 

.35 

5 

.18 

.11 

.23 

.42 

1.00 

.71 

.43 

.59 

6 

-.17 

.38 

.44 

.86 

.71 

1.00 

.45 

.53 

7 

.14 

.38 

.21 

.57 

.43 

.45 

1.00 

.72 

8 

.01 

-.03 

.29 

.35 

.59 

.53 

.72 

1.00 

Recall  that  r   is  the  estimated  serial  correlation  for  a 
n 

simulated  Normal  (0,1)  subsample  of  size  n.   Also,  the 

estimated  correlations  shown  above  were  computed  using  10 

2  2 

values  (replications)  of  s  (r   )  and  s  (r   )  for  each  i  and 

n .  n . 

i  1 


44 


all  i  and  j ,  where  the  estimation  was  done  by  repeating  the 
SIMTB1  experiment  with  10  different  batches  of  50,000  simu- 
lated random  variables.   Since  only  10  values  went  into  each 
correlation  calculation,  the  table  is  only  accurate  to 
within  approximately  ±  2//T0  =  .632.   We  see  some  indication 
of  positive  correlation,  especially  when  i  and  j  are  close, 
but  the  lack  of  independence  is  not  severe  enough  to  hurt 
the  regression  results  for  either  the  estimated  means  or 
variances  significantly. 

A  second  assumption,  implicit  in  any  regression,  is  that 
the  dependent  variables  have  equal  variances.   This  condition 
holds  true  for  the  means,  which  can  be  shown  to  satisfy 


Var(0(ni))   =   | 


independently  of  i.   The  estimated  variances,  however,  are  not 
equivalent  and,  if  we  assume  the  0.(n.)  to  be  approximately 


Normally  distributed  so  that 


I  N/n.    . 
L£   y  (0.  (n.)  -  0(n.)) 
j-1      J      X       X 


2 
is  approximately  proportional  to  a  xM  hi/n   i   i  random  variable, 

L  V 

we  can  compute 


Var(s2(n. 


MNn . -n2 
l   l 
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2  . 

To  correct  this  problem,  SIMTB1  scales  the  s  (n.)  by  /nT 

so  that 


. 2  2      ~   2 

Var ( /n .  s  ( n . ) )   = 


l      i  '      MN  -  n .      MN 

l 

since  MN  >>  n..  The  design  matrix  is  scaled  accordingly  and 
the  values  b_  ,  b,  ,  b2,  and  b^  discussed  in  Section  VLB.  are 
reported. 

Table  2  shows  the  effects  of  the  rescaling  by  presenting 

2 
first  the  estimated  variances  of  the  s  (n.) ,  where  the  esti- 

i 

mation  is  done  by  repeating  SIMTB1  for  10  batches  of  50,000 
simulated  data  points.   These  estimated  variances  decrease 
as  n.  increases,  closely  paralleling  the  second  line  of  Table 

2  which  has  the  approximate  theoretical  values  (i.e., 

2 
2/(MNn.  -n.))  .   The  final  line  of  Table  2  shows  the  estimated 
l    l 

2  2 

variances  of  the  rescaled  s  (n.)#  i.e.,  the  /n .  s  (n.), 

l  i     i 

which,  as  expected  and  hoped,  show  a  more  constant  variance 
with  i. 

Although  future  versions  of  SIMTB1  will  include  more 
sophisticated  regression  routines  and  the  ability  to  gener- 
ate independent  samples  at  each  subsample  size,  the  SIMTB1 
is  quick,  usable,  and  accurate  for  most  situations. 

D.   INTERPRETING  THE  SERIAL  CORRELATION  RESULTS 

Returning  to  Figure  3(a)  which  shows  the  simulated  dis- 
tribution of  the  serial  correlation  coefficient  from  indepen- 
dent, Normal  (0,1)  data,  the  following  comments  summarize  the 
most  striking  features: 


TABLE  2 

A  COMPARISON  OF  THE  ESTIMATED  VARIANCE  OF  s2(r   ) 

WITH  THE  APPROXIMATE  THEORETICAL  VARIANCE  OF   ni 

s  (r   )  AND  WITH  THE  APPROXIMATELY  EQUIVARIANT  SCALED 

i        VERSIONS,  n~-5  s2(r   ). 

l        n. 

1         5 
All  entries  have  been  multiplied  by  10  . 


n.  =  10    20    30    40    50    75   100   150 

l 

Var(s2(r  ))  .177  .150  .204  .079   .047  .031  .049  .022 

n. 

l 


Approx.  Theoretical 
"n 


Var(s2(r  ))  .400  .200  .133  .100  .080  .053  .040  .027 


l 

Var(nT*5  s2(r  ))        1.77  2.99  6.12  3.18  2.33  2.33  4.88  3.35 
i      n. 

l 


2  2 

The  estimated  variances  of  s  (r   )  and  /n.  s  (r   )  were 

n .         i     n . 

1  X2 

calculated  using  10  independent  replications  of  s  (r   ) . 

i 
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(a)  The  boxplots  appear  very  symmetric  at  all  subsample 
sizes  with  nearly  equal  numbers  of  outliers  at  either  tail 
and  with  mean  and  median  coincidental.   This  observation  is 
confirmed  by  the  estimates  of  skewness  in  the  summary  sta- 
tistics.  Kurtosis  is  mildly  negative  at  small  subsample 
sizes  but,  overall,  asymptotic  Normality  seems  to  take  hold 
rather  quickly. 

(b)  The  average  serial  correlation  is  negative  for  small 
subsam-les.   This  is  demonstrated  by  the  dotted  regression 
curve  which  starts  at  approximately  -.10  and  levels  off  near 

0  for  subsamples  greater  than  about  85.   The  dashed  asymptote 
of  -.000272  is  very  close  to  the  theoretical  value  of  0,  and 
the  mean  values  in  the  summary  table  closely  reflect  the  bias 
of  -1/n. 

(c)  The  standard  deviations  of  the  simulated  distributions 
are  very  close  to  the  asymptotic  values  of  n.  *  ,  although 
the  lead  coefficient  in  the  regression  on  the  variances, 

bQ  =  .743756,  is  not  as  close  to  the  theoretical  value  of  1 
as  we  would  hope.   When  SIMTB1  is  repeated  10  times  with  10 
different  batches  of  simualted  data,  we  find  an  average  value 
for  bn  to  be  1.0604,  with  a  standard  deviation  for  bQ  of 
.307.   The  estimation  procedure  for  bQ ,  therefore,  remains 
valid,  but  the  estimate  itself  is  highly  variable. 

The  agreement  between  the  simulated  and  the  theoretical, 
asymptotic  values  of  the  bias  and  variances  was  discovered 
previously  by  Cox  (19  66) .   SIMTB1  has  now  allowed  us  to 
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automatically  look  at  a  broader  range  of  subsample  sizes  and 
to  see,  through  boxplots  and  estima  es  of  skewness  and 
kurtosis,  a  fuller  picture  of  any  changes  in  the  distribution 
of  the  estimator.   We  can  be  satisfied  that  estimates  of 
serial  correlation  do  behave  approximately  as  Normal  (-1/n, 
1/n)  rangome  variables  when  the  underlying  data  is  Normal 
(0,1)  . 

If  the  lead  terms  in  the  expansions  of  the  mean  and 
variance  of  the  estimated  correlation  coefficient  (i.e.,  a. , 
a,,  and  b-.)  had  been  unknown,  we  would  also  have  a  fairly 

good  idea  now  of  what  they  were. 

2 
When  the  underlying  data  is  Xn t    Figure  3(b)  confirms 

Cox's  observation  that  the  bias  is  relatively  uneffected  but, 

for  small  subsamples,  the  standard  deviation  is  smaller  than 

-1/2 
the  expected  n   '  .   Unlike  Figure  3(a) ,  there  is  a  pronounced 

skewness  in  the  boxplots  in  Figure  3(b)  with  many  more  out- 
liers at  the  positive  end,  and  with  the  mean  higher  than  the 
median  at  the  first  four  subsample  sizes.   The  problem  of 
suppressed  variance  seems  cured  at  n_  =  100  and  ng  =  150, 
but  the  skewness  remains  and  could  cause  problems  in  tests 
of  significance. 

Figure  3(c),  which  is  based  on  an  underlying  batch  of 
simulated  Lognormal  (0,1)  data,  shows  a  slight  exaggeration 
of  the  effects  in  Figure  3(b).   The  standard  deviation  is 
more  suppressed  and  does  not  attain  the  theoretical  level 
by  nQ  =  150.   The  positive  skewness  is  more  pronounced  and 

o 

kurtosis  does  not  approach  the  theoretical  value  of  0. 
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Overall,  the  effects  of  long-tailed  data  on  the  distribu- 
tion of  the  serial  correlation  coefficient  can  be  summarized 
as  follows: 

(i)     Bias  is  not  significantly  effected  and  remains  at 

approximately  -1/n. 
(ii)    The  variance  of  the  distribution  of  the  serial 
correlation  coefficient  is  reduced  by  longer- 
tailed  data. 
(iii)   Positive  skewness  is  created  in  the  distribution, 
(iv)    Kurtosis  may  become  positive  at  larger  subsample 

sizes  . 
(v)     For  long-tailed  data  (i.e.,  Lognormal) ,  a  subsample 
size  of  150  is  not  large  enough  to  insure 
asymptotic  Normality. 

E.   SIMTB1  OUTPUT  FOR  THE  Z-TRANSFORM  OF  THE  CORRELATION 

Fisher's  z-transform  of  the  estimated  correlated  coeffi- 
cient is  defined  by: 


i      1  +  r 

Z    =   7T  log  t , 

n     2    3  1  -  r 

n 


where  r   is  the  estimated  serial  correlation  presented  in 
n  r 

Section  VLB.   The  transformation  is  intended  to  make  the 

distribution  of  the  Z   more  Normal  than  that  of  the  r  . 

n  n 

When  the  same  SIMTBl  experiment  described  in  Section  VLB. 
is  run  using  Z  as  the  estimator  instead  of  r  ,  we  get  the 
results  shown  in  Figures  4(a),  4(b)  and  4(c).   It  should 
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be  noted  that  the  scale  of  the  boxplots  here  has  been  forced 

to  be  approximately  comparable  to  the  scale  for  the  boxplots 

in  Section  VI. B.   This  is  done  by  suppressing  outliers  that 

are  more  than  1.5  interquartile  distances  beyond  the  quar- 

tiles  of  the  boxplot.   If  we  had  allowed  the  data  to  scale 

the  boxplots,  we  would  have  seen  a  much  wider  range  on  the 

vertical  axis  because  the  Z   are  not  restricted  to  the  limits 

n 

of  -1  to  +1  and  because  there  is  one  far  outlier  at  -3.8. 
In  this  type  of  "reduced  graphics ,"  we  still  see  the  number 
of  outliers  that  fall  beyond  the  allowable  range  through  the 
numbers  at  the  ends  of  the  boxplots,  but  we  do  not  see  their 
actual  locations. 

Figure  4 (a)  shows  the  distribution  of  the  z-transformed 
correlation  coefficients  when  the  underlying  data  is  simu- 
lated, Normal  (0,1) .   At  each  subsample  size,  the  mean  and 

standard  deviation  are  close  to  the  theoretical  n   and 

-1/2 
n     respectively.   The  skewness  and  kurtosis  at  subsample 

size  n,  =  10  are  far  from  the  theoretical  Normal  distribution 

values  of  0  and  0,  reflecting  partly  the  one  far  outlier  at 

-3.8  and  partly  the  negative  skew  in  the  remainder  of  the 

Z   's.   For  other  subsample  sizes,  there  is  no  strong  evi- 

nl 
dence  to  contradict  the  assumption  of  approximate  Normality. 

The  relationship  between  Figure  4(b)  and  4(a)  is  similar 

to  that  between  3(b)  and  3(a).   Figure  4(b),  which  is  based 

on  simulated  Xi  data,  shows  (a)  bias  that  is  the  same  as  for 

the  transformed  correlations  based  on  Normal  data,  (b)  slightly 
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suppressed  variances,  particularly  at  small  subsample  sizes 
and  (c)  positive  skewness  which  persists  at  large  subsample 
sizes.   In  addition,  there  are  signs  of  positive  kurtosis 
at  small  subsample  sizes. 

Figure  4(c)  is  based  on  Lognormal  (0,1)  data  and  shows 
high  values  of  skewness  and  kurtosis  at  almost  all  subsample 
sizes.   Approximate  Normality  seems  an  unwarranted  assumption. 
In  fact,  the  kurtosis  is  converging  very  slowly  to  its 
asymptotic  value  of  0 . 

In  general,  using  the  z-transform  does  not  help  with 
Normality  assumptions,  especially  when  dealing  with  long- 
tailed  distributions. 

F.   SIMTB1  OUTPUT  FOR  THE  2-FOLD  JACKKNIFE  OF  THE  CORRELATION 

The  final  Figures,  5(a),  5(b)  and  5(c),  deal  with  the  2- 
fold  jackknife  estimate  of  correlation.   Again,  the  figures 
are  reduced  graphics  with  scaling  comparable  to  that  of  the 
boxplots  of  Sections  VI. D  and  VI. E.   To  define  the  estimator, 
we  start  with  a  given  subsample  of  size  n,  compute  the  serial 
correlation  for  the  first  I n/2J  points  and  call  it  r, (n/2), 
compute  the  serial  correlation  for  the  second  jn/2|  points  and 
call  it  r2(n/2)  and  compute  the  serial  correlation  for  the 
entire  subsample  of  n  points  and  call  it  rQ (n) .   Each  compu- 
tation follows  the  formula  in  Section  VI. b.   The  three  esti- 
mators are  then  combined  to  form  two  pseudo-values, 


r1*(n)      =      2rQtn)    -    r±Ln/2) 
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and 


r2*(n)   =   2rQ(n)  -  r2(n/2)  , 


and  the  final  jackknife  estimator  for  that  subsample  is  de- 
fined as 


r,*n)  +  r2*(n) 
r(n)   =   = 


Although  a  jackknife  estimator  may  have  many  favorable 
properties,  we  are  concerned  here  primarily  with  its  ability 
to  remove  bias,  hopefully  without  inflating  the  variances 
of  the  estimator  and/or  inducing  nonnormality . 

Figure  3(a) ,  based  on  simulated  Normal  (0,1)  data,  shows 
nearly  complete  removal  of  bias,  even  at  small  subsample 
sizes.  .  The  cost  of  the  bias  reduction  is  reflected  in  an  in- 
crease of  nearly  50%  in  the  standard  deviation  of  the  corre- 
lation estimate  for  subsample  size  10,  and  lesser  relative 
increases  at  larger  subsample  sizes.   There  is  also  an  indica- 
tion of  a  positive  skew  for  small  subsample  sizes,  and  the 
problem  that  the  jackknife  estimator  need  not  fall  into  the 
-1  to  +1  range  which  is  desirable  for  a  correlation  coeffi- 
cient estimate. 

2 
When  using  simulated  x-i  data  as  in  Figure  5(b),  or  simu- 
lated Lognormal  (0,1)  data  as  in  Figure  5(c),  there  is  again 
no  problem  with  bias.   Variance  inflation,  though  it  exists 
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at  small  subsample  sizes,  is  not  as  large  as  when  Normal 
(0,1)  data  is  used.   The  distributions  of  the  jackknifed 
correlations  show  very  pronounced  positive  skews,  however,  as 
well  as  positive  kurtosis.   These  two  problems  are  worse  for 
the  longer-tailed  Lognormal  data. 

Overall,  the  jackknife  estimator  is  very  successful  at 
removing  bias  but  the  costs  include  variance  inflation,  which 
can  be  severe  at  small  subsample  sizes,  plus  increased  posi- 
tive skewness  and  kurtosis  when  the  estimates  are  based  on 
data  from  longer-tailed  distributions. 

G.   COMPARISON  OF  THE  THREE  ESTIMATES  OF  CORRELATION 

For  Normal  (0,1)  data,  the  distribution  of  the  usual 
correlation  coefficient  displayed  in  Figure  3(a)  behaves 
very  much  as  theoretical  asymptotic  calculations  would  pre- 
dict, even  at  small  subsample  sizes.   This  makes  it  possible 
to  correct  for  bias  in  the  estimator  and  to  perform  tests 
of  significance.   Use  of  Fisher's  z-transform,  as  illustrated 
in  Figure  4(a)  does  not  seem  necessary  since  it  does  not 
significantly  improve  the  approximate  Normality  of  the  esti- 
mator.  The  jackknife  estimator  in  Figure  5(a)  may  be  valuable 
if  a  direct,  unbiased  estimator  is  needed  but  the  inflated 
variance  of  the  jackknife  estimator  may  limit  the  usefulness 
of  the  estimate  as  well  as  make  any  tests  of  significance 
too  conservative. 

When  the  underlying  data  comes  from  a  longer-tailed  dis- 
tribution, the  usual  correlation  coefficient  in  Figures  3(b) 
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and  3(c)  retains  a  predictable  bias  term,  although  the  vari- 
ance of  its  distribution  is  slightly  depressed  and  the  skew- 
ness  and  kurtosis  becomes  positive,  even  for  subsamples  as 
large  as  150.   This  means  that  it  is  still  possible  to 
estimate  the  correlation  accurately,  but  tests  of  signifi- 
cance fall  on  shakey  assumptions  of  Normality.   The  z-transform 
in  Figures  4(b)  and  4(c)  does  little  to  firm  up  those  assump- 
tions and,  in  some  cases,  makes  the  situation  worse.   As  in 
the  case  of  Normal  data,  the  2-fold  jackknifed  correlation  in 
Figures  5(b)  and  5(c)  is  bias-free  but  follows  a  fairly  non- 
Normal  distribution  which  would  invalidate  significance 
testing. 

All  of  the  preceding  observations  and  conclusions  flow 
immediately  from  the  nine  figures  presented  so  far.   Further 
studies  could  easily  be  done  through  SIMTB1,  looking  at  larger 
subsample  sizes,  correlated  data,  and  alternative  marginal 
distributions.   For  demonstration  purposes,  though,  it  is 
better  to  proceed  to  our  second  application. 
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VII.   STUDY  OF  PROBLEMS  OF  ESTIMATING  SHAPE 

PARAMETERS  FOR  HIGHLY  SKEWED  DISTRIBUTIONS 

A.   ESTIMATING  THE  SHAPE  PARAMETER  FOR  A  GAMMA  DISTRIBUTION 

As  a  second  application  of  SIMTB1,  we  will  consider  a 
problem  which  has  received  much  less  statistical  attention; 
asymptotic  results  are  summarized  in  Cox  and  Lewis  (19  66, 
Ch.  2)  [Ref.  7]  and  Johnson  and  Kotz  (1970,  Ch.  17)  [Ref. 
8].   We  want  to  estimate  the  shape  parameter,  K,  for  a  Gamma 
distribution,  where  the  Gamma  density  is  given  by 


.  K.  K  x      -Kx/u        .  n    „  .  n      .  n 
( — )   -fTT — r  e  x>0;   K>0;   y  >  0 

u    1  (x 


f  (x)   = 


x  <  0 


Notice  that  the  mean  of  this  distribution  is  u,  not  K/u  as 
in  some  differently  parameterized  versions  of  the  Gamma  den- 
sity.  For  the  data  that  will  be  simulated  for  use  in  SIMTB1 
we  will  use  K  =  5  and  u  =  1  and  K  =  0.25  and  jj  =  1.   The 
closer  the  mean  of  our  estimate  is  to  5  or  0.25,  the  better 
(in  terms  of  bias)  is  our  estimation  procedure.   Other  factors 
such  as  the  variance  and  Normality  of  the  estimator  will  of 
course  also  have  influence  in  the  determination  of  a  preferred 
estimator. 

Section  VII. B  will  compare  the  commonly  used  maximum 
likelihood  estimator  to  the  competing  method  of  moments 
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estimator.   Both  procedures  result  in  asymptotically  Normal 
estimators  (Cramer,  1948)  but  the  m.l.e.  is  usually  pre- 
ferred because  of  its  favorable  asymptotic  relative  efficiency 
(Cox  and  Lewis,  1966)  [Ref.  7].   Through  SIMTBl,  though,  we 
will  see  that  for  small  subsamples  the  estimated  variances 
of  the  two  estimators  of  K  are  not  as  far  apart  as  asymp- 
totic results  lead  us  to  believe.   In  addition,  the  bias  that 
appears  in  both  estimators  is  smaller  for  the  moment 
estimator . 

In  Section  VII. C.  we  will  use  a  four-fold  jackknife  of 
both  the  m.l.e.  and  moments  estimators  to  successfully  remove 
the  bias.   What  is  remarkable  is  that,  unlike  the  jackknifing 
of  the  serial  correlation,  there  is  little  or  no  cost  in 
terms  of  variance  inflation  and  nonnormality  for  the  jack- 
knifed  moment  estimator.   When  K  =  .25,  we  will  see  in  Section 
VII. D.  that  the  jackknifed  m.l.e.  dominates  the  other  three 
estimators  at  all  subsample  sizes  when  considering  the  mean, 
variance,  and  Normality  of  the  estimator. 

B.   MAXIMUM  LIKELIHOOD  AND  MOMENT  ESTIMATORS  OF  K 

Figure  6(a)  is  very  similar  in  format  to  the  figures 
that  have  already  been  presented  for  the  correlation  example 
except  that: 

CI)   The  estimator  whose  distribution  is  being  displayed 
is  the  maximum  likelihood  estimator  of  K,  the  shape  parameter 
of  a  Gamma (5)  population.   We  denote  the  estimator,  computed 
from  a  simulated  subsample  of  size  n,  by  K(n)  and  define  it 
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to  be  the  solution  of  the  equation; 


/v  n  n 

n[log  K(n)  -V(K(n))]   =   n  log   £   X./n  -      log  x.  , 

i=l   x     i=l      x 


where  the  X.  are  the  simulated  Gamma (5)  random  variables 
and  ¥(•)  is  the  digamma  function  (Cox  and  Lewis,  1966). 

(2)  The  eight  subsample  sizes  which  we  will  be  looking  at 

are  n,  =  33,  n.  =  50,  n_  =  71,  n.  =  100,  nc  =  125,  nr    =  166, 
1  2.  6  4  5  o 

n_  =  250  and  ng  =  500.   We  will  not  see  as  much  detail  at 
small  subsample  sizes  but  we  will  see  some  of  the  asymptotic 
(n  =  500)  effects  coming  in. 

(3)  At  each  subsample  size  we  will  work  with  M*  =  20 
independent  replications  of  N*  =  2500  simulated  Gamma (5) 
random  variables,  instead  of  the  M  =  10  replications  of 
N  =  5000  variables  used  previously.   The  total  number  of 
independent  simulated  random  variables  across  replications 
remain  constant  at  the  program  maximum  of  5  0,000.   Hence, 
the  boxplot  at  subsample  size  50  in  Figure  6(a)  represents 
the  distribution  of  M*  [n*/50J  ■  1000  estimates  of  K(50)  just 
as  the  boxplot  at  subsample  size  50  in  Figure  3 (a)  represents 
the  distribution  of  M  |n/50|  =  1000  estimates  of  r(n).   As 
long  as  the  product,  Mxn,  remains  constant,  the  only  effect 
that  changing  the  number  of  replications  has,  up  to  rounding 
in  |N/n.|  ,  is  to  change  the  results  in  the  regression  on  the 
averages.   By  using  M*  =  20  and  N*  =  2500,  SIMTB1  reports 
regression  coefficients  averaged  over  20  replications,  but, 
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within  each  replication,  the  dependent  variables  are  averages 
over  just  12  50  0/n .  I  values  of  the  estimator. 

(4)   The  boxplots  are  presented  using  the  reduced  graphics 
option.   In  this  option  any  extreme  outliers  (i.e.,  those 
beyond  1.5  interquartile  distances)  are  included  as  a  count 
at  the  tail  of  each  boxplot.   This  option  was  chosen  in  order 
to  give  more  graphical  weight  to  the  body  of  the  distributions 
and  the  fall-off  in  the  bias.   Limited  printer  resolution 
makes  it  impossible  to  show  details  in  the  body  and  the  tails 
of  the  distributions  if  there  are  many  straggling  outliers. 
In  the  case  of  very  extreme  outliers,  no  detail  would  be  seen 
in  the  body  of  the  boxplot  without  the  reduced  graphics  option, 

Figure  6(b)  looks  at  the  distribution  of  the  moment 
estimator  of  K,  the  shape  parameter  of  a  Gamma  (K)  population: 

-?   n       -  2 
K(n)   =   (n-1)  XV  I       (X.  -X)   , 

i=l    x 

n 
where   X  =   T   X./n.  n  is  the  subsample  size,  and  the  X. 

i=l   x 
are  the  simulated  Gamma (5)  random  variables.   The  SIMTB1 

options  and  parameters  mentioned  in  (2),  (3)  and  (4)  preceding 

are  also  in  effect  here. 

The  two  Figures,  6(a)  and  6(b),  show  a  very  pronounced 

bias  in  both  estimation  procedures,  although  the  moment 

estimator  is  slightly  closer  to  the  unbiased  value  of  5. 

As  expected,  the  standard  deviation  of  the  m.l.e.  is  lower 

than  that  of  the  moment  estimator  although  the  relative 
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difference  at  small  subsample  sizes,  for  instance  1.448  versus 
1.482  at  n,  =  33 ,  may  not  outweight  the  increase  in  bias  with 
the  m.l.e.   At  larger  subsample  sizes,  the  relative  differ- 
ence is  close  to  the  theoretical  asymptotic  relative  efficiency 
of  .78  (i.e.,  .91  at  n   =  250). 

Both  estimators  also  show  distributions  with  positive 
skewness  and  kurtosis  that  decrease  to  the  asymptotic  0 
levels  as  subsample  size  increases.   The  asymptotics  appear 
to  take  hold  more  quickly  for  the  moment  estimator  than  for 
the  m.l.e. 

In  summary,  SIMTB1  shows  that  the  m.l.e.  is  indeed  better 
than  the  moment  estimator  in  terms  of  variance,  but  not  as 
good  for  small  sample  sizes  as  asymptotic  results  would  lead 
us  to  believe.   In  the  other  areas  of  bias  and  asymptotic 
Normality,  the  moment  estimator  would  have  to  be  preferred. 

C.   4-FOLD  JACKKNIFED  ESTIMATORS  OF  K 

Figures  6(c)  and  6(d)  show  the  distributions  of  the  4-fold 
jackknife  m.l.e.  of  K  and  4-fold  jackknifed  moment  estimator 
of  K,  respectively.   A  4-fold  jackknife  estimator  is  similar 
to  the  2-fold  jackknife  estimator  described  in  Section  VI. F. 
except  that  there  are  4  pseudo-values  that  come  out  of  dividing 
each  subsample  into  fourths.   More  details  can  be  found  in 
Mosteller  and  Tukey  (1977)  [Ref .  6] . 

The  purpose  of  the  jackknife  is  to  remove  the  conspicuous 
bias  observed  in  Figures  6(a)  and  6(b).  This  goal  is  seen  to 
be  accomplished  in  Figure  6(c)  and  6(d)  and  we  can  also  note 
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smaller  values  of  skewness  and  kurtosis,  indicating  a  quicker 
approach  to  asymptotic  Normality.   The  skewness  and  kurtosis 
of  the  jackknifed  moment  estimator  are  the  lowest,  at  small 
subsample  sizes,  among  all  estimators.   The  variance  of  the 
jackknifed  moment  estimator  is  also  only  slightly  inflated, 
as  is  the  variance  of  the  jackknifed  m.l.e. 

All  told,  the  jackknifed  moment  estimator,  because  of 
its  lack  of  bias,  small  variance,  and  low  skewness  and  kurtosis, 
would  be  the  method  of  choice  if  estimation  of  K  or  signifi- 
cance testing  was  the  goal. 

D.  RESULTS  FOR  K  =  0.2  5 

In  Figures  7(a),  7(b),  7(c)  and  7(d)  we  show  similar 
results  to  those  discussed  above  for  the  case  K  =  5.0,  but 
using  K  =  0.25.   The  fact  (Cox  and  Lewis,  19  66,  Ch.  3)  [Ref. 
7]  that  the  m.l.e.  estimate  is  much  more  efficient  than  the 
moment  estimate  is  graphically  illustrated.   What  is  new  is 
the  effect  of  jackknifing:   bias  is  reduced  without  the 
sacrifice  of  variance  inflation  or  nonnormality . 

Further  comparisons  and  interpretations  are  similar  to 
those  done  for  the  case  K  =  5.0,  and  are  left  to  the  reader. 

E.  CONCLUSIONS 

Simply  by  providing  SIMTB1  with  the  desired  estimators, 
we  have  been  able  (a)  to  explore  in  depth  the  effects  of 
changes  in  data  distribution  and  of  different  estimation 
procedures  on  the  calculation  of  the  serial  correlation 
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coefficient,  and  (b)  to  compare  four  different  ways  to  esti- 
mate the  shape  parameter  in  a  highly  skewed  Gamma  population 

The  graphics  and  numerical  output  combine  to  let  us  see 
and  quantify  distributional  changes  that  occur  as  subsample 
size  grows.   We  can  see  bias  fall  away,  variance  shrink,  and 
skewness  disappear  as  the  estimator  approaches  asymptotic 
Normality.   Terms  in  the  asymptotic  expansion  of  the  mean 
and  variance  of  the  estimator  are  automatically  calculated 
and  can  be  used  to  compare  different  estimators. 
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VIII.   COMPARISON  OF  DIFFERENT  METHODS  FOR  ESTIMATING 
THE  VARIABILITY  OF  THE  STANDARD  DEVIATION  OF 
CORRELATION  ESTIMATES 

A.  INTRODUCTION 

Bradley  Efron  and  Gail  Gong  (1982)  review  in  their 
article,  "A  Leisurely  Look  at  the  Bootstrap,  the  Jackknife, 
and  Cross-validation"  [Ref.  9]  different  methods  for  esti- 
mating statistical  error  of  parameter  estimates.   They  dis- 
cuss in  particular  the  problem  of  estimating  the  error  of  a 
statistical  estimator  with  the  example  of  the  estimates  of 
the  standard  error  for  the  correlation  coefficient  of  a 
bivariate  normal  distribution.   They  compared  the  standard 
deviation  estimates  of  different  methods  using  200  simulations 
at  one  fixed  sample  size  of  14. 

SIMTB2  was  used  to  explore  the  distributions  of  the 
estimates  they  used  in  their  article.   The  bootstrap,  the 
jackknife,  the  infinitesimal  jackknife  (delta  method)  and 
the  normal  theory  are  the  methods  for  which  the  distributions 
were  explored.   The  estimation  methods  will  not  be  explained 
in  detail  (see  Efron  and  Gong  in  The  American  Statistician, 
Feb.  1983,  Vol.  37,  No.  1),  only  the  setup  of  SIMTB2  and 
the  output  will  be  discussed. 

B.  SETUP  OF  SIMTB2  AND  THE  ESTIMATOR  FUNCTIONS 

The  bivariate  normal  distributed  input  data  with  known 
correlation  0.5  was  generated  with  an  IMSL  random  number 
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generator  (GGNSM) .   14,000  bivariate  data  points  were  gener- 
ated and  then  sectioned  into  10  blocks,  each  having  1,400 
points  (M  =  10,  N  =  1400)."   The  subsample  sizes  (NE(k))  used 
are  10,  14,  20,  28,  35,  40,  70  and  100.   Only  subsample  size 
14  was  used  by  Efron  and  Gong. 

The  program  had  to  be  run  with  5  different  standard 
deviation  estimator  functions.   SIMTB2,  as  the  other  ver- 
sions, can  only  handle  up  to  3  estimator  functions  in  one 
program  run.   To  make  the  given  outputs  comparable  the  fixed 
scale  option  (SVS  =  1,  YMIN  =0.0  and  YMAX  =  0.7)  was  chosen. 

The  bootstrap  estimate  of  standard  deviation  was  done  with 
2  different  numbers  of  bootstrap  replications  (B  =  128, 
B  =  512) .   So  for  the  subsample  size  of  14,  the  bootstrap 
procedure  with  B  bootstrap  replications  in  itself  was  done 
(1400/14  -  100;  100*10  =  1000)  1,000  times. 

The  jackknife  function  followed  the  standard  jackknife 
procedure  and  used  the  jackknife  formula  for  the  standard 
deviation: 


n  ~  1/2 

n-1   r   ,tt      -    ,2 


jack       n   i^1    (l)     (  •) 


with 


0(j)   =   G(xrx2 xi-l'xi+l'""V   for   X(i) 


and 
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V)  ■  5  X   °<i)   for  x(.) 


The  delta  method  (infinitesimal  jackknife)  function 
calculated  the  estimate  with  the  formula: 


Delta  '  l4nL-2     -2"-      -2      "   *    "  ~   ~   J ' 
y20    y02    y20y02    ylx     ylly02    ^11^02 

with 


U9h 


=   |  I     (Y  -Y)g(Z  -Z)h 

n  i=l   x       x 


and 


X.   =   (Y.,Z.) 


In  doing  the  actual  calculations  for  the  delta  method, 
for  small  subsample  sizes  ( 10  ,14  ,  . . . , 40)  the  value  of  the 
variance  becomes  negative.   A  negative  variance  can  not  be 
interpreted  with  much  meaning.   The  number  of  negative  values 
goes  down  (from  305  for  NE(1)  =  10  to  18  for  NE(6)  =  40) 
with  increasing  subsample  size.   Most  of  the  negative  values 
are  small.   To  solve  the  programming  problem  (square  root 
of  a  negative  number)  the  function  sets  this  negative  value 
to  0.0D0. 

For  the  normal  theory  estimate,  instead  of  carrying  out 
a  bootstrap,  an  approximation  formula  was  used.   It  is  the 
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same  formula  the  article  uses  for  the  comparison 


~2 

S       =   (1  ~    P  > 
norm 


^(n  -  3) 
See  Efron  and  Gong  for  details. 

C.   INTERPRETING  THE  SIMTB2  OUTPUT 

The  output  of  the  program  runs  is  provided  as  Figures  6, 
7,  8,  9  and  10.   A  comparison  of  the  numerical  output  (sub- 
sample  size  14)  with  the  results  is  done  in  Table  3. 

The  bootstrap  procedure  (Figures  8  and  9)  was  done  with 
2  different  numbers  of  bootstrap  replications  (B  =  128 
and  B  =  512) .   Both  distributions  for  the  standard  deviation 
(S.D.)  look  very  similar.   Both  are  positively  skewed  with 
some  outliers  at  the  right  tail.   In  both  cases  the  outliers 
are  in  the  same  range.   As  the  boxplots  and  the  summary 
statistics  show,  the  increase  of  the  number  of  bootstrap 
replications  (B)  does  not  result  in  a  large  improvement  in 
the  performance  of  the  estimation  function. 

The  jackknife  estimate  (Figure  10)  has  a  positively  skewed 
distribution  with  outliers.   The  distribution  of  the  jack- 
knife  estimate  looks  very  similar  to  the  distribution  of  the 
bootstrap  estimates.   For  small  subsample  sizes  the  bootstrap 
distribution  has  more  outliers.   Overall  the  performance  of 
the  jackknife  procedure  is  as  good  as  the  bootstrap,  but  the 
jackknife  needs  less  computer  time. 
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The  distribution  of  the  estimates  produced  by  the  delta 
method  (Figure  11)  is  negatively  skewed  and  has  nearly  no 
outliers.   But  in  calculating  the  estimates  the  problem  of 
negative  values  for  the  variance  came  up.   For  some  sub- 
samples,  the  final  estimate  (the  standard  deviation)  could 
not  be  calculated,  since  the  corresponding  value  of  the 
variance  was  negative.   In  these  cases,  the  standard  devia- 
tion was  set  to  0.0.   This  procedure  influences  the  distribu- 
tion and  the  summary  statistics.   The  influence  is  more 
important  for  small  subsample  sizes  than  for  larger  ones. 
So  the  graphical  and  numerical  output  should  be  seen  with 
this  fact  in  mind. 

The  normal  theory  function  (Figure  12)  produces  esti- 
mates with  a  negatively  skewed  distribution  but  only  a  few 
outliers  and  the  distribution  is  tailed  to  the  left.   The 
tail  of  the  distribution  is  in  the  opposite  direction  of 
all  other  distributions.   For  the  estimate  of  the  standard 
deviation  for  the  correlation  coefficient  the  result  of  the 
normal  theory  is  close  to  the  true  result.   This  may  not  be 
valid  for  other  estimators. 

In  addition  to  the  comparisons  Efron  and  Gong  did,  with 
SIMTB2  it  is  easy  to  investigate  how  the  sample  size  will 
influence  the  estimate  of  the  standard  deviation  of  the  corre- 
lation coefficient.   In  Table  4  the  methods  are  compared  for 
a  subsample  size  of  10  and  100.   With  increasing  subsample 
size  the   quality  of  the  estimate  should  improve,  but  the 
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improvement  may  be  different  for  the  different  methods  of 
estimation.   By  making  the  subsample  size  10  times  larger, 
with  the  SIMTB2-side  effect  of  reducing  the  total  number  of 
calculated  estimates,  the  bootstrap  improves  less  than  the 
jackknife,  delta  method  and  normal  theory. 
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IX.   FINAL  CONCLUSIONS 

SIMTBED  with  the  different  versions  can  be  used  on 
digital  computers  of  different  size  (mainframe  to  micro) 
and  type.   The  limitations  in  using  the  program  are  given 
more  by  hardware  constraints  like  memory  size  and  computer 
time  than  by  the  program  itself. 

The  FORTRAN  program  is  completely  portable,  changes  in 
the  code,  may  only  be  necessary  to  adapt  the  program  to 
special  restrictions  given  by  a  special  type  of  hardware. 
This  may  occur  in  using  micro  computers  more  often  than 
with  mainframe  computers.   Up  to  now  all  versions  of  code 
are  written  for  the  more  normal  standard  computer  environment 
and  do  not  need  special  equipment  (color  printer,  etc.) .  * 
Additionally  hardware  dependent  features  like  color  output 
can  improve  the  graphics  of  the  program. 

SIMTBED  makes  it  easy  to  evaluate  the  result  of  statis- 
tical experiments.   The  combination  of  graphics  and  numerical 
summaries  for  different  sample  sizes  make  it  easy  to  judge 
the  distributional  behavior  of  a  statistical  estimator.   The 
result  can  be  seen  without  additional  computations  in  the 
graphic  outputs.   Comparing  only  the  baxplots  it  is  possible 
to  judge  the  influence  of  subsample  size  on  the  variability 
of  an  estimator. 
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Besides  for  research  the  program  can  be  used  in  showing 
students  the  distributional  behavior  of  different  estimators 
in  a  pictorial  way.   It  is  easy  to  compare  the  different 
behavior  of  similar  estimators  (e.g.,  biased  vs.  unbiased 
estimator  of  the  variance)  for  different  sample  sizes. 

The  easy  use  and  the  fast  visual  impression  of  the 
distributional  behavior  of  an  estimator/  given  by  the 
graphic  output  is  one  of  the  advantages  in  using  SIMTBED. 
Besides  this  fast  first  visual  impression  all  necessary  and 
needed  numerics  are  given  for  further  and  deeper 
investigations . 
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